{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 쇼어(Shor)의 소인수 분해 알고리즘 써큐 예제"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이 노트북은 쇼어의 알고리즘을 써큐에서 구현하는 교육 자료입니다. 이 지침서는 [이 곳의 써큐 예제](https://github.com/quantumlib/Cirq/blob/master/examples/shor.py)에서 수정/보완된 자료입니다.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"노트북을 위한 라이브러리 가져오기.\"\"\"\n",
    "import fractions\n",
    "import math\n",
    "import random\n",
    "\n",
    "import numpy as np\n",
    "import sympy\n",
    "from typing import Callable, List, Optional, Sequence, Union\n",
    "\n",
    "import cirq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 위수 찾기"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "정수 $n$을 소인수 분해하는 문제는 <i>모듈러 지수함수</i>(차후에 설명)의 주기를 찾는 문제로 좁혀질 수 있습니다. 이 주기를 찾는 것은 (매우 높은 확률로) 법(modulo) $n$에 대한 곱셈 군(multiplicative group)의 무작위로 추출된 한 원소의 <i>위수</i>를 찾는 것으로 달성할 수 있습니다.\n",
    "\n",
    "양의 정수 $n$에 대하여, \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbb{Z}_n := \\{x \\in \\mathbb{Z}_+ : x < n \\text{ 이고 } \\text{gcd}(x, n) = 1\\}\n",
    "\\end{equation}\n",
    "\n",
    "인 법 $n$에 대한 곱셈 군을 정의합니다.\n",
    "주어진 $x \\in \\mathbb{Z}_n$에 대하여 $x^r \\text{ mod } n = 1$을 만족하는 가장 작은 양의 정수 $r$를 계산합니다. \n",
    "\n",
    "군이론/정수론으로부터 다음을 보일 수 있습니다.\n",
    "\n",
    "(1) 그러한 정수 $r$이 존재합니다. (Note that $g^{|G|} = 1_G$ for any group $G$ with cardinality $|G|$ and element $g \\in G$., but it's possible that $r < |G|$.)\n",
    "\n",
    "(2) 만일 소수 $p$와 $q$에 대하여 $n = pq$일 때, $|\\mathbb{Z}_n| = \\phi(n) = (p - 1) (q - 1)$입니다. (이러한 함수 $\\phi$를 [오일러의 피 함수(Euler's totient function)](https://ko.wikipedia.org/wiki/%EC%98%A4%EC%9D%BC%EB%9F%AC_%ED%94%BC_%ED%95%A8%EC%88%98)라고 합니다.)\n",
    "\n",
    "(3) 모듈러 지수함수\n",
    "\n",
    "\\begin{equation}\n",
    "f_x(z) := x^z \\mod n\n",
    "\\end{equation}\n",
    "\n",
    "는 주기 $r$ (원소 $x \\in \\mathbb{Z}_n$의 위수)을 가집니다. 즉, $f_x(z + r) = f_x(z)$입니다. \n",
    "\n",
    "(4) 모듈러 지수함수의 주기를 알 수 있다면, (매우 높은 확률로) $p$와 $q$를 알아낼 수 있습니다. -- 즉, $n$의 인수분해입니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "복습차원에서, 정수 $n$의 곱셈군 $\\mathbb{Z}_n$의 원소들을 다음의 단순한 함수로 그려볼 수 있습니다.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"곱셈군 Z_n의 원소들을 계산하는 함수.\"\"\"\n",
    "def multiplicative_group(n: int) -> List[int]:\n",
    "    \"\"\"법 n의 곱셈군을 반환합니다.\n",
    "    \n",
    "    Args:\n",
    "        n: 곱셈군의 법(modulus).\n",
    "    \"\"\"\n",
    "    assert n > 2\n",
    "    group = [1, 2]\n",
    "    for x in range(3, n):\n",
    "        if math.gcd(x, n) == 1:\n",
    "            group.append(x)\n",
    "    return group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "예를 들어, 법 $n = 15$의 곱셈군은 아래와 같습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The multiplicative group modulo n = 15 is:\n",
      "[1, 2, 4, 7, 8, 11, 13, 14]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"곱셈군 예제.\"\"\"\n",
    "n = 15\n",
    "print(f\"The multiplicative group modulo n = {n} is:\")\n",
    "print(multiplicative_group(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "평범한 곱셈 연산에서 이러한 원소들의 집합이 정말로 군을 형성하는지 확인해 보세요."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 고전적 위수 찾기"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "아래에 정수 $x \\in \\mathbb{Z}_n$의 위수 $r$을 고전적으로 찾는 함수가 있습니다. 이 함수는 단순히 수열 \n",
    "\n",
    "\\begin{align}\n",
    "    &x^2 \\text{ mod } n, \\\\\n",
    "    &x^3 \\text{ mod } n, \\\\\n",
    "    &x^4 \\text{ mod } n, \\\\\n",
    "    &\\ \\ \\ \\ \\ \\ \\ \\ \\vdots\n",
    "\\end{align}\n",
    "\n",
    "을 $x^r = 1 \\text{ mod } n$인 정수 $r$이 될 때까지 계산합니다. $|\\mathbb{Z}_n| = \\phi(n)$이기 때문에, 이 위수 찾기 알고리즘은 시간복잡도 $O(\\phi(n))$을 갖는데 이는 비효율적입니다. (정수 $n$의 비트수를 $L$이라 할 때 대략 $O(2^{L / 2})$.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Z_n의 한 원소의 위수를 고전적으로 계산하는 함수.\"\"\"\n",
    "def classical_order_finder(x: int, n: int) -> Optional[int]:\n",
    "    \"\"\"x**r mod n == 1를 만족하는 최소의 양의 정수를 계산.\n",
    "\n",
    "    Args:\n",
    "        x: 위수가 계산될 1보다 크고 법 n의 곱셈군에 속하는 정수. (정수 n과 서로소인 양의 정수\n",
    "           들의 곱으로 이루어진 수)\n",
    "        n: 상기 곱셈군의 법\n",
    "\n",
    "    Returns:\n",
    "        x**r == 1 mod n를 만족하는 최소의 양의 정수.\n",
    "        이 알고리즘은 항상 성공합니다. (즉, 절대로 None을 반환하지 않습니다.)\n",
    "\n",
    "    Raises:\n",
    "        ValueError: x가 1이거나 법 n의 곱셈군의 원소가 아닐 때.\n",
    "    \"\"\"\n",
    "    # x가 Z_n에 속하는 유효한 값인지 확인합니다.\n",
    "    if x < 2 or x >= n or math.gcd(x, n) > 1:\n",
    "        raise ValueError(f\"Invalid x={x} for modulus n={n}.\")\n",
    "    \n",
    "    # 위수를 결정합니다.\n",
    "    r, y = 1, x\n",
    "    while y != 1:\n",
    "        y = (x * y) % n\n",
    "        r += 1\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "주어진 $x \\in \\mathbb{Z}_n$와 $n$에 대하여 위수 $r$을 계산하는 예제가 아래 코드 영역에 있습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x^r mod n = 8^4 mod 15 = 1\n"
     ]
    }
   ],
   "source": [
    "\"\"\"(고전적으로) 한 원소의 위수를 계산하는 예제.\"\"\"\n",
    "n = 15  # 곱셈군은 [1, 2, 4, 7, 8, 11, 13, 14]\n",
    "x = 8\n",
    "r = classical_order_finder(x, n)\n",
    "\n",
    "# 위수가 정말로 맞는지 확인합니다.\n",
    "print(f\"x^r mod n = {x}^{r} mod {n} = {x**r % n}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "하지만 쇼어 알고리즘의 핵심 양자 요소는 바로 위수 찾기, 양자 회로로 구현되어야 하기에 아래에서 논의하도록 하겠습니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 양자 위수 찾기"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "양자 위수 찾기는 본질적으로 무작위로 추출된 $x \\in \\mathbb{Z}_n$에 대한 모듈러 지수함수 $f_x(z)$를 계산하는 유니타리 $U$를 수반하는 양자 위상 추정 알고리즘 입니다. $U$가 기초 게이트들로부터 어떻게 계산되는지의 상세한 이론은 특히 처음 읽는 경우 설명하기 복잡할 수 있습니다. 이 강의 자료에서는, 기초 게이트들에 대한 상세한 이론을 들여다 보지 않고도 그러한 유니타리 $U$ 구현할 수 있는 써큐의 산술연산을 사용할 것입니다.\n",
    "\n",
    "아래는 써큐에서 단순 산술 연산(덧셈)의 예제를 처음 보여줍니다. 그후 우리가 다룰 연산(모듈러 지수)을 논의할 것입니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 써큐에서 양자 산술 연산하기"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "여기서 우리는 써큐에서 산술 연산, 즉 모듈러 덧셈(modular addition)을 정의하는 예제에 대해 논의할 것입니다. 이 연산은 입력 레지서터의 값을 목적 레지스터에 더한 뒤 나머지를 취해 넣습니다. 더 자세히 설명하면 이 연산은 두 큐비트 레지스터에서 다음과 같이 동작합니다.\n",
    "\n",
    "\\begin{equation}\n",
    "|a\\rangle_i |b\\rangle_t \\mapsto |a\\rangle_i |a + b \\text{ mod } N_t \\rangle_t .\n",
    "\\end{equation}\n",
    "\n",
    "여기서 첨자 $i$와  $t$는 각각 영어 <i>i</i>nput(입력)과 <i>t</i>arget(목적) 레지스터를 의미합니다. 그리고 $N_t$는 목적 레지스터의 차원을 의미합니다.\n",
    "\n",
    "이 연산을 정의하기 위해 (이름은 `Adder`라고 합시다) `cirq.ArithmeticOperation` 클래스에서 상속받아 아래의 4개 메소드 함수를 재정의 합니다. 주요 함수인 `apply`는 산술을 정의합니다. 여기서 우리는 위에 나온 더 정확한 $a + b \\text{ mod } N_t$ 대신에 단순히 수식 $a + b$을 구현해 봅시다. -- `cirq.ArithmeticOperation` 클래스는 그 연산이 반드시 가역적이어야 하기 때문에 우리가 단순히 $a + b$로 무엇을 의미하려 했는지 유추할 수 있습니다. . "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"써큐에서 (양자) 산술 연산자 정의하는 예제.\"\"\"\n",
    "class Adder(cirq.ArithmeticOperation):\n",
    "    \"\"\"양자 덧셈.\"\"\"\n",
    "    def __init__(self, target_register, input_register):\n",
    "        self.input_register = input_register\n",
    "        self.target_register = target_register\n",
    "    \n",
    "    def registers(self):\n",
    "        return self.target_register, self.input_register\n",
    "    \n",
    "    def with_registers(self, *new_registers):\n",
    "        return Add(*new_registers)\n",
    "    \n",
    "    def apply(self, target_value, input_value):\n",
    "        return target_value + input_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이제 클래스를 정의했으므로 회로에서 사용할 수 있습니다. 아래는 두 개 큐비트 레지스터를 생성하고 \n",
    "$X$게이트를 이용해 첫 번째 레지스터를 (이진수) $|10\\rangle$로, 두 번째 레지스터를 (이진수)\n",
    "$|01\\rangle$로 설정합니다. . \n",
    "그런 후 `Adder`연산을 사용하고, 모든 큐비트를 측정합니다.\n",
    "\n",
    "이진수로 $10 + 01 = 11$이기 때문에, 목적 레지스터에서 항상 $|11\\rangle$를 측정할 것이라 기대할 수 있습니다.\n",
    "게다가, 입력 레지스터를 바꾸지 않았기 때문에 입력 레지스터도 항상 $|10\\rangle$를 측정할 것입니다. \n",
    "요약하면, 우리가 측정할 수 있는 것은 $1011$ 뿐임을 알 수 있습니다. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Circuit:\n",
      "\n",
      "0: ───X───#3──────────────────────────────────────────M───\n",
      "          │\n",
      "1: ───────#4──────────────────────────────────────────M───\n",
      "          │\n",
      "2: ───────<__main__.Adder object at 0x7ff2b8c3d9b0>───M───\n",
      "          │\n",
      "3: ───X───#2──────────────────────────────────────────M───\n",
      "\n",
      "\n",
      "Measurement outcomes:\n",
      "\n",
      "   0  1  2  3\n",
      "0  1  0  1  1\n",
      "1  1  0  1  1\n",
      "2  1  0  1  1\n",
      "3  1  0  1  1\n",
      "4  1  0  1  1\n"
     ]
    }
   ],
   "source": [
    "\"\"\"회로에서 Adder를 사용하는 예제.\"\"\"\n",
    "# 두 개의 큐비트 레지스터\n",
    "qreg1 = cirq.LineQubit.range(2)\n",
    "qreg2 = cirq.LineQubit.range(2, 4)\n",
    "\n",
    "# 회로 정의하기\n",
    "circ = cirq.Circuit(\n",
    "    cirq.ops.X.on(qreg1[0]),\n",
    "    cirq.ops.X.on(qreg2[1]),\n",
    "    Adder(input_register=qreg1, target_register=qreg2),\n",
    "    cirq.measure_each(*qreg1),\n",
    "    cirq.measure_each(*qreg2)\n",
    ")\n",
    "\n",
    "# 화면에 회로 출력하기\n",
    "print(\"Circuit:\\n\")\n",
    "print(circ)\n",
    "\n",
    "# 측정 결과 출력하기\n",
    "print(\"\\n\\nMeasurement outcomes:\\n\")\n",
    "print(cirq.sample(circ, repetitions=5).data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이 코드 영역의 출력을 살펴보면 가장 먼저 초기화용 $X$ 게이트들을 볼 수 있고, `Adder` 연산 그리고 마지막 측정들을 볼 수 있습니다. 다음으로 우리가 예측한대로 모든 측정 비트열들이 $1011$임을 알 수 있습니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "덧셈 연산의 유니타리를 보는것도 아래 처럼 가능합니다. 여기서 우리는 목적 레지스터를 \n",
    "$|00\\rangle$를 갖는 2개 큐비트로 설정했습니다. 입력 레지스터는 $|01\\rangle$에 해당하도록\n",
    "1개 큐비트를 1로 설정하였습니다. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ryan/programs/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:5: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  \"\"\"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, 1],\n",
       "       [1, 0, 0, 0],\n",
       "       [0, 1, 0, 0],\n",
       "       [0, 0, 1, 0]], dtype=int32)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Adder 연산의 유니타리 계산 예제.\"\"\"\n",
    "cirq.unitary(\n",
    "    Adder(target_register=cirq.LineQubit.range(2),\n",
    "          input_register=1)\n",
    ").astype(np.int32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이 유니타리를 다음과 같이 이해할 수 있습니다. 유니타리의 $i$번째 열은 상태 $|i + 1 \\text{ mod } 4\\rangle$를 의미합니다. \n",
    "예를들어 $0$번째 열을 보면 $|i + 1 \\text{ mod } 4\\rangle = |0 + 1 \\text{ mod } 4\\rangle = |1\\rangle$임을 알 수 있습니다. \n",
    "똑같이 $1$번째 열을 보면, $|i + 1 \\text{ mod } 4\\rangle = |1 + 1 \\text{ mod } 4\\rangle = |2\\rangle$ 상태입니다. \n",
    "마찬가지로 나머지 두 개 열도 이해할 수 있습니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 모듈러 지수 산술 연산."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "모듈러 지수 산술 연산도 아래 나와 있는 것처럼 단순 덧셈 산술 연산과 비슷한 방법으로 정의할 수 있습니다. \n",
    "쇼어 알고리즘을 이해하기 위한 목적으로, 다음 코드 영역에서 가장 중요한 요소는 산술 연산을 정의하는 `apply` 메소드 입니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"쇼어 알고리즘에서 사용되는 모듈러 지수 연산을 정의하기.\"\"\"\n",
    "class ModularExp(cirq.ArithmeticOperation):\n",
    "    \"\"\"양자 모듈러 거듭제곱.\n",
    "\n",
    "    이 클래스는 밑을 지수 만큼 곱한 뒤 주어진 법에 따라 모듈로 연산을 하는 유니타리를 표현합니다. \n",
    "    더 정확하게는, 모듈러 지수 x**e mod n를 계산하는 유니타리 V를 표현합니다.\n",
    "\n",
    "        V|y⟩|e⟩ = |y * x**e mod n⟩ |e⟩     0 <= y < n\n",
    "        V|y⟩|e⟩ = |y⟩ |e⟩                  n <= y\n",
    "\n",
    "    여기서 y는 목적 레지스터이고, e는 지수 레지스터, x는 밑, n은 법입니다. 결론적으로,\n",
    "\n",
    "        V|y⟩|e⟩ = (U**e|r⟩)|e⟩\n",
    "\n",
    "    인 유니타리 U는 다음과 같이 정의됩니다.\n",
    "\n",
    "        U|y⟩ = |y * x mod n⟩      0 <= y < n\n",
    "        U|y⟩ = |y⟩                n <= y\n",
    "    \"\"\"\n",
    "    def __init__(\n",
    "        self, \n",
    "        target: Sequence[cirq.Qid],\n",
    "        exponent: Union[int, Sequence[cirq.Qid]], \n",
    "        base: int,\n",
    "        modulus: int\n",
    "    ) -> None:\n",
    "        if len(target) < modulus.bit_length():\n",
    "            raise ValueError(f'Register with {len(target)} qubits is too small '\n",
    "                             f'for modulus {modulus}')\n",
    "        self.target = target\n",
    "        self.exponent = exponent\n",
    "        self.base = base\n",
    "        self.modulus = modulus\n",
    "\n",
    "    def registers(self) -> Sequence[Union[int, Sequence[cirq.Qid]]]:\n",
    "        return self.target, self.exponent, self.base, self.modulus\n",
    "\n",
    "    def with_registers(\n",
    "            self,\n",
    "            *new_registers: Union[int, Sequence['cirq.Qid']],\n",
    "    ) -> cirq.ArithmeticOperation:\n",
    "        if len(new_registers) != 4:\n",
    "            raise ValueError(f'Expected 4 registers (target, exponent, base, '\n",
    "                             f'modulus), but got {len(new_registers)}')\n",
    "        target, exponent, base, modulus = new_registers\n",
    "        if not isinstance(target, Sequence):\n",
    "            raise ValueError(\n",
    "                f'Target must be a qubit register, got {type(target)}')\n",
    "        if not isinstance(base, int):\n",
    "            raise ValueError(\n",
    "                f'Base must be a classical constant, got {type(base)}')\n",
    "        if not isinstance(modulus, int):\n",
    "            raise ValueError(\n",
    "                f'Modulus must be a classical constant, got {type(modulus)}')\n",
    "        return ModularExp(target, exponent, base, modulus)\n",
    "\n",
    "    def apply(self, *register_values: int) -> int:\n",
    "        assert len(register_values) == 4\n",
    "        target, exponent, base, modulus = register_values\n",
    "        if target >= modulus:\n",
    "            return target\n",
    "        return (target * base**exponent) % modulus\n",
    "\n",
    "    def _circuit_diagram_info_(\n",
    "            self,\n",
    "            args: cirq.CircuitDiagramInfoArgs,\n",
    "    ) -> cirq.CircuitDiagramInfo:\n",
    "        assert args.known_qubits is not None\n",
    "        wire_symbols: List[str] = []\n",
    "        t, e = 0, 0\n",
    "        for qubit in args.known_qubits:\n",
    "            if qubit in self.target:\n",
    "                if t == 0:\n",
    "                    if isinstance(self.exponent, Sequence):\n",
    "                        e_str = 'e'\n",
    "                    else:\n",
    "                        e_str = str(self.exponent)\n",
    "                    wire_symbols.append(\n",
    "                        f'ModularExp(t*{self.base}**{e_str} % {self.modulus})')\n",
    "                else:\n",
    "                    wire_symbols.append('t' + str(t))\n",
    "                t += 1\n",
    "            if isinstance(self.exponent, Sequence) and qubit in self.exponent:\n",
    "                wire_symbols.append('e' + str(e))\n",
    "                e += 1\n",
    "        return cirq.CircuitDiagramInfo(wire_symbols=tuple(wire_symbols))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`apply` 메소드에서 `(target * base**exponent) % modulus`를 계산하고 있습니다. \n",
    "`target`과 `exponent` 변수는 각각의 큐비트 레지스터의 값들에 의존합니다.\n",
    " 그리고 밑 `base`와 법 `modulus`는 상수로, `modulus`는 $n$이고 `base`는 $x \\in \\mathbb{Z}_n$입니다. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "우리가 사용할 총 큐비트 개수는 $3 (L + 1)$인데 $L$은 인수분해할 정수 $n$을 저장하기 위한 총 비트열의 길이입니다. \n",
    "그러므로 모듈러 지수를 계산할 유니타리의 크기는 $4^{3(L + 1)}$입니다. \n",
    "적당히 크지 않은 수 $n = 15$에 대해서도 유니타리는 $2^{30}$개의 부동소수점 실수들을 메모리에 저장해야합니다.\n",
    "이는 대부분의 표준 노트북 컴퓨터에서는 엄두내지 못할 크기입니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "To factor n = 15 which has L = 4 bits, we need 3L + 3 = 15 qubits.\n"
     ]
    }
   ],
   "source": [
    "\"\"\"위상 추정을 위해 목적레지스터와 지수 레지스터를 생성하고\n",
    "쇼어 알고리즘에서 필요한 큐비트 개수를 살펴봅니다.\n",
    "\"\"\"\n",
    "n = 15\n",
    "L = n.bit_length()\n",
    "\n",
    "# 목적 레지스터는 L개의 큐비트를 갖습니다.\n",
    "target = cirq.LineQubit.range(L)\n",
    "\n",
    "# 지수 레지스터는 2L + 3개 큐비트를 갖습니다.\n",
    "exponent = cirq.LineQubit.range(L, 3 * L + 3)\n",
    "\n",
    "# 정수 n을 인수분해 하기 위한 총 큐비트 개수를 화면에 출력합니다.\n",
    "print(f\"To factor n = {n} which has L = {L} bits, we need 3L + 3 = {3 * L + 3} qubits.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "단순 덧셈 연산과 비교하면, 이 모듈러 지수 연산은 다음과 같이 (메모리가 허용하는)출력할 수 있는 유나타리를 갖습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"모듈러 지수 연산을 위한 유니타리(의 일부분) 확인하기.\"\"\"\n",
    "# 법 n의 곱셈군의 한 원소를 고릅니다.\n",
    "x = 5\n",
    "\n",
    "# 유니타리의 일부분을 확인합니다. n이 충분히 작은 경우에만 주석을 해제합시다.\n",
    "# cirq.unitary(ModularExp(target, exponent, x, n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 회로에서 모듈러 지수 연산을 사용하기"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "쇼어 알고리즘의 양자 요소는 모듈러 지수 연산에 대응하는 유니타리 $U$의 위상추정입니다.\n",
    "다음 코드 영역에서는 위에서 정의한 `ModularExp`를 사용하는 쇼어 알고리즘 회로를 생성합니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"위수 찾기 양자 회로를 생성하는 함수.\"\"\"\n",
    "def make_order_finding_circuit(x: int, n: int) -> cirq.Circuit:\n",
    "    \"\"\"x 모듈로 n의 위수를 계산하는 양자 회로를 반환합니다.\n",
    "\n",
    "    이 회로는 양자 위상 추정을 사용하여 다음 유니타리의 고유값을 계산합니다.\n",
    "\n",
    "        U|y⟩ = |y * x mod n⟩      0 <= y < n\n",
    "        U|y⟩ = |y⟩                n <= y\n",
    "\n",
    "    Args:\n",
    "        x: 법 n에 대해 위수를 찾을 대상인 밑이 되는 양의 정수\n",
    "        n: x의 위수와 서로소인 법\n",
    "\n",
    "    Returns:\n",
    "        x 모듈로 n의 위수를 찾는 양자 알고리즘 회로.\n",
    "    \"\"\"\n",
    "    L = n.bit_length()\n",
    "    target = cirq.LineQubit.range(L)\n",
    "    exponent = cirq.LineQubit.range(L, 3 * L + 3)\n",
    "    return cirq.Circuit(\n",
    "        cirq.X(target[L - 1]),\n",
    "        cirq.H.on_each(*exponent),\n",
    "        ModularExp(target, exponent, x, n),\n",
    "        cirq.QFT(*exponent, inverse=True),\n",
    "        cirq.measure(*exponent, key='exponent'),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이 함수를 사용하여, 주어진 $x$와 $n$에 대한 회로를 그려볼 수 있습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: ────────ModularExp(t*7**e % 15)────────────────────────────\n",
      "           │\n",
      "1: ────────t1─────────────────────────────────────────────────\n",
      "           │\n",
      "2: ────────t2─────────────────────────────────────────────────\n",
      "           │\n",
      "3: ────X───t3─────────────────────────────────────────────────\n",
      "           │\n",
      "4: ────H───e0────────────────────────QFT^-1───M('exponent')───\n",
      "           │                         │        │\n",
      "5: ────H───e1────────────────────────#2───────M───────────────\n",
      "           │                         │        │\n",
      "6: ────H───e2────────────────────────#3───────M───────────────\n",
      "           │                         │        │\n",
      "7: ────H───e3────────────────────────#4───────M───────────────\n",
      "           │                         │        │\n",
      "8: ────H───e4────────────────────────#5───────M───────────────\n",
      "           │                         │        │\n",
      "9: ────H───e5────────────────────────#6───────M───────────────\n",
      "           │                         │        │\n",
      "10: ───H───e6────────────────────────#7───────M───────────────\n",
      "           │                         │        │\n",
      "11: ───H───e7────────────────────────#8───────M───────────────\n",
      "           │                         │        │\n",
      "12: ───H───e8────────────────────────#9───────M───────────────\n",
      "           │                         │        │\n",
      "13: ───H───e9────────────────────────#10──────M───────────────\n",
      "           │                         │        │\n",
      "14: ───H───e10───────────────────────#11──────M───────────────\n"
     ]
    }
   ],
   "source": [
    "\"\"\"주기 찾기 양자 회로의 예제.\"\"\"\n",
    "n = 15\n",
    "x = 7\n",
    "circuit = make_order_finding_circuit(x, n)\n",
    "print(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이전에 설명했듯이, 아다마르 게이트를 통해 지수 레지스터를 동등한 중첩상태를 만듭니다.\n",
    "목적레지스터의 마지막 큐비트에 가한 $X$는 위상 반동에 사용됩니다.\n",
    "모듈러 지수 연산은 위상추정에서 제어 유나타리들의 나열들을 수행합니다.\n",
    "그 후, 역 양자 푸리에 변환을 지수 레지스터에 가한 뒤 결과를 측정해 읽어내면 됩니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "측정 결과를 설명하기 위해, 더 작은 회로에서 추출을 해보겠습니다. \n",
    "(실제로는 절대로 쇼어 알고리즘을 짝수인 $n = 6$에서 수행하지 않습니다. 이 것은 단순히 측정 결과를 설명하기 위한 예제일 뿐 입니다.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Raw measurements:\n",
      "exponent=01010001, 00000000, 00000000, 00000000, 00000000, 00000000, 00000000, 00000000, 00000000\n",
      "\n",
      "Integer in exponent register:\n",
      "   exponent\n",
      "0         0\n",
      "1       256\n",
      "2         0\n",
      "3       256\n",
      "4         0\n",
      "5         0\n",
      "6         0\n",
      "7       256\n"
     ]
    }
   ],
   "source": [
    "\"\"\"쇼어의 주기 찾기 회로를 측정하기.\"\"\"\n",
    "circuit = make_order_finding_circuit(x=5, n=6)\n",
    "res = cirq.sample(circuit, repetitions=8)\n",
    "\n",
    "print(\"Raw measurements:\")\n",
    "print(res)\n",
    "\n",
    "print(\"\\nInteger in exponent register:\")\n",
    "print(res.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "측정된 각각의 비트열들을 정수로 해석할 수 있습니다. 하지만 이 정수들이 우리에게 무엇을 말해줄까요?\n",
    "다음 소단원에서는 이들을 해석하기위해 고전적으로 어떻게 후처리를 하는지 볼 것입니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 고전적 후처리"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "우리가 측정한 정수는 $x \\in \\mathbb{Z}_n$의 위수 $r$에 대해 $s / r$에 가깝습니다. \n",
    "이때 $0 \\le s < r$인 정수 입니다. $s / r$로 부터 $r$을 구하기 위해 \n",
    "연분수 알고리즘(continued fractions algorithm)을 사용합니다. \n",
    "만일 위수 찾기 회로가 성공한다면 그 값을 출력하고, 아니면 `None`값을 출력합니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_measurement(result: cirq.TrialResult, x: int, n: int) -> Optional[int]:\n",
    "    \"\"\"위수 찾기 회로의 출력을 해석하기.\n",
    "\n",
    "    구체적으로, exp(2πis/r)이 다음 유니타리 U의 고유값인 s/r을 결정합니다.\n",
    "\n",
    "        U|y⟩ = |xy mod n⟩  0 <= y < n\n",
    "        U|y⟩ = |y⟩         n <= y\n",
    "    \n",
    "    그 후 가능하다면 (연분수 알고리즘으로) r을 계산하고 반환합니다.\n",
    "\n",
    "    Args:\n",
    "        result: make_order_finding_circuit에 의해 생성된 회로의 출력을 추출하여 얻은\n",
    "                시행 결과.\n",
    "\n",
    "    Returns:\n",
    "        r(x 모듈로 n의 위수) 혹은 None.\n",
    "    \"\"\"\n",
    "    # 지수 레지스터의 출력 정수를 읽습니다.\n",
    "    exponent_as_integer = result.data[\"exponent\"][0]\n",
    "    exponent_num_bits = result.measurements[\"exponent\"].shape[1]\n",
    "    eigenphase = float(exponent_as_integer / 2**exponent_num_bits)\n",
    "\n",
    "    # f = s / r를 결정하기 위한 연분수 알고리즘.\n",
    "    f = fractions.Fraction.from_float(eigenphase).limit_denominator(n)\n",
    "    \n",
    "    # 분자가 0이면 위수 찾기가 실패 했으므로 None을 반환합니다.\n",
    "    if f.numerator == 0:\n",
    "        return None\n",
    "    \n",
    "    # 찾은 분모가 실제 위수이면 그 값을 반환합니다.\n",
    "    r = f.denominator\n",
    "    if x**r % n != 1:\n",
    "        return None\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "다음 코드 영역은 위수 찾기 회로를 생성하기 실행한 뒤 고전적 후처리를 통해 위수를 계산하는\n",
    "예제를 나타냅니다. 이 알고리즘의 양자 요소는 확률적으로 성공함을 상기해봅시다.\n",
    " 만일 위수가 `None`이라면 몇번 더 코드영역을 재실행 해보세요."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finding the order of x = 5 modulo n = 6\n",
      "\n",
      "Raw measurements:\n",
      "exponent=1, 0, 0, 0, 0, 0, 0, 0, 0\n",
      "\n",
      "Integer in exponent register:\n",
      "   exponent\n",
      "0       256\n",
      "\n",
      "Order r = 2\n",
      "x^r mod n = 5^2 mod 6 = 1\n"
     ]
    }
   ],
   "source": [
    "\"\"\"고전적 후처리 예제.\"\"\"\n",
    "# n과 x를 설정합니다.\n",
    "n = 6\n",
    "x = 5\n",
    "\n",
    "print(f\"Finding the order of x = {x} modulo n = {n}\\n\")\n",
    "measurement = cirq.sample(circuit, repetitions=1)\n",
    "print(\"Raw measurements:\")\n",
    "print(measurement)\n",
    "\n",
    "print(\"\\nInteger in exponent register:\")\n",
    "print(measurement.data)\n",
    "\n",
    "r = process_measurement(measurement, x, n)\n",
    "print(\"\\nOrder r =\", r)\n",
    "if r is not None:\n",
    "    print(f\"x^r mod n = {x}^{r} mod {n} = {x**r % n}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "곱셈군 $\\mathbb{Z}_6$의 원소 $x = 5$의 위수가 $r = 2$임을 알 수 있습니다.\n",
    " 실제로, $5^2 \\text{ mod } 6 = 25 \\text{ mod } 6 = 1$로 확인 가능합니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 양자 위수 측정기"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이제 우리가 지금까지 작성한 함수들을 사용하여 위수 찾기의 양자 버전을 이용한 효율적인 함수를 정의할 수 있습니다. \n",
    "아래의 양자 위수 측정기는 회로를 생성하고 실행한 뒤 측정 결과를 처리합니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantum_order_finder(x: int, n: int) -> Optional[int]:\n",
    "    \"\"\"x**r mod n == 1인 최소의 양의 정수 r을 계산한다.\n",
    "    \n",
    "    Args:\n",
    "        x: 위수가 계산될 정수. 항상 1보더 크며 법 n의 곱셈군에 속한다. 이때 x는 n과 서로소인\n",
    "           양의 정수들로 이루어져 있습니다.\n",
    "        n: 위 곱셈군의 법\n",
    "    \"\"\"\n",
    "    # 정수 x가 법 n 곱셈군의 유효한 원소인지 확인합니다.\n",
    "    if x < 2 or n <= x or math.gcd(x, n) > 1:\n",
    "        raise ValueError(f'Invalid x={x} for modulus n={n}.')\n",
    "\n",
    "    # 위수 찾기 회로를 생성합니다.\n",
    "    circuit = make_order_finding_circuit(x, n)\n",
    "    \n",
    "    # 위수 찾기 회로에서 결과를 추출합니다.\n",
    "    measurement = cirq.sample(circuit)\n",
    "    \n",
    "    # 측정결과를 처리한 후 반환합니다.\n",
    "    return process_measurement(measurement, x, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "이제 우리가 구현한 위수 측정기의 양자 구현과 쇼어 알고리즘의 양자 요소를 완료합니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 완성된 인수 분해 알고리즘"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "양자 위수 측정기(또는 고전적 위수 측정기)를 사용해 쇼어의 알고리즘을 완성할 수 있습니다.\n",
    " 다음 코드 영역에서는 몇가지 전처리 단계를 추가합니다.\n",
    "\n",
    "(1) $n$이 짝수인지 확인하기\n",
    "\n",
    "(2) $n$이 소수인지 확인하기\n",
    "\n",
    "(3) $n$이 소수의 거듭제곱인지 확인하기\n",
    "\n",
    "위 세가지 모두 고전 컴퓨터에서 효율적으로 처리 가능합니다. 덧붙여 \n",
    " we add the last necessary post-processing step which uses the order  \n",
    " 마지막으로 위수 $r$을 사용하여 $n$의 자명하지 않은 인수인 $p$를 계산하는 후처리 단계를 추가합니다. \n",
    " 이는 ($r$이 짝수라고 가정하면) $y = x^{r / 2} \\text{ mod } n$을 계산하여 얻어낼 수 있습니다.\n",
    "  그렇다면 $p = \\text{gcd}(y - 1, n)$입니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" 인수분해의 시작부터 끝까지하는 함수.\"\"\"\n",
    "def find_factor_of_prime_power(n: int) -> Optional[int]:\n",
    "    \"\"\"n이 소수 한개의 거듭제곱인 경우 자명하지 않은 n의 인수를 반환합니다.\n",
    "    그 외에는 None을 반환합니다.\"\"\"\n",
    "    for k in range(2, math.floor(math.log2(n)) + 1):\n",
    "        c = math.pow(n, 1 / k)\n",
    "        c1 = math.floor(c)\n",
    "        if c1**k == n:\n",
    "            return c1\n",
    "        c2 = math.ceil(c)\n",
    "        if c2**k == n:\n",
    "            return c2\n",
    "    return None\n",
    "\n",
    "\n",
    "def find_factor(\n",
    "    n: int,\n",
    "    order_finder: Callable[[int, int], Optional[int]] = quantum_order_finder,\n",
    "    max_attempts: int = 30\n",
    ") -> Optional[int]:\n",
    "    \"\"\"합성수 n의 자명하지 않은 인수를 반환합니다.\n",
    "\n",
    "    Args:\n",
    "        n: 인수분해할 정수.\n",
    "        order_finder: 법 n 곱셈군에 속하는 원소의 위수를 찾는 함수.\n",
    "        max_attempts: 함수 order_finder가 호출되는 횟수의 상한값이자 시도해볼 무작위수\n",
    "                      x의 개수.\n",
    "\n",
    "    Returns:\n",
    "        n의 자명하지 않은 인수 혹은 그러한 인수가 없을 때 None.\n",
    "        n의 인수 k는 1이거나 n일때 자명합니다.\n",
    "    \"\"\"\n",
    "    # 만일 n이 소수이면 자명하지 않은 인수는 없습니다.\n",
    "    if sympy.isprime(n):\n",
    "        print(\"n is prime!\")\n",
    "        return None\n",
    "    \n",
    "    # 짝수이면 자명하지 않은 인수는 2입니다.\n",
    "    if n % 2 == 0:\n",
    "        return 2\n",
    "    \n",
    "    # n이 소수의 거듭제곱이라면, 효율적으로 자명하지 않은 인수를 구할 수 있습니다.\n",
    "    c = find_factor_of_prime_power(n)\n",
    "    if c is not None:\n",
    "        return c\n",
    "    \n",
    "    for _ in range(max_attempts):\n",
    "        # 무작위로 2에서 n - 1 사이의 값을 고릅니다.\n",
    "        x = random.randint(2, n - 1)\n",
    "        \n",
    "        # x와 n이 서로소 인지 확인합니다.\n",
    "        c = math.gcd(x, n)\n",
    "        \n",
    "        # x와 n이 서로소가 아니라면, 운이 좋게도 한 번에 비자명한 인자를 구해냈습니다.\n",
    "        if 1 < c < n:\n",
    "            return c\n",
    "        \n",
    "        # 법 n으로 하는 x의 위수 r을 위수 측정기로 계산합니다.\n",
    "        r = order_finder(x, n)\n",
    "        \n",
    "        # 위수 측정기가 실패하면 재시도합니다.\n",
    "        if r is None:\n",
    "            continue\n",
    "        \n",
    "        # 위수가 홀수이면 재시도합니다.\n",
    "        if r % 2 != 0:\n",
    "            continue\n",
    "        \n",
    "        # 비자명한 인자를 계산합니다.\n",
    "        y = x**(r // 2) % n\n",
    "        assert 1 < y < n\n",
    "        c = math.gcd(y - 1, n)\n",
    "        if 1 < c < n:\n",
    "            return c\n",
    "\n",
    "    print(f\"Failed to find a non-trivial factor in {max_attempts} attempts.\")\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "함수 `find_factor`는 `quantum_order_finder`를 기본으로 사용하며 쇼어의 알고리즘을 사용하게 됩니다. \n",
    "이전에 언급했듯이, 고전적으로 이 회로를 시뮬레이션하는데에는 매우 큰 메모리가 필요하기 때문에\n",
    "$n \\ge 15$인 정수에 대해서는 쇼어의 알고리즘을 실행할 수 없습니다. 이때는 대체제로 고전적인 위수\n",
    "찾기 알고리즘을 사용할 수 있습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Factoring n = pq = 184573\n",
      "p = 487\n",
      "q = 379\n"
     ]
    }
   ],
   "source": [
    "\"\"\"쇼어 알고리즘(위수 측정기)을 통해 인수분해하기 예제.\"\"\"\n",
    "# 인수분해할 수를 정의합니다.\n",
    "n = 184573\n",
    "\n",
    "# 인수를 찾아봅니다.\n",
    "p = find_factor(n, order_finder=classical_order_finder)\n",
    "q = n // p\n",
    "\n",
    "print(\"Factoring n = pq =\", n)\n",
    "print(\"p =\", p)\n",
    "print(\"q =\", q)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}